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Abstract 



We investigate a population of primordial binaries during the ini- 
tial stage of evolution of a star cluster. For our calculations we assume 
that equal mass stars form rapidly in a tidally truncated gas cloud, 
that ~ 10% of the stars are in binaries and that the resulting star 
cluster undergoes an epoch of violent relaxation. We study the col- 
lisional interaction of the binaries and single stars, in particular, the 
ionization of the binaries and the energy exchange between binaries 
and single stars. We find that for large N systems (iV > 10'^) even the 
most violent beginning leaves the binary distribution function largely 
intact. Hence, the binding energy originally tied up in the cloud's 
protostellar pairs is preserved during the relaxation process and the 
binaries are available to interact at later times within the virialized 
cluster. 

Subject headings: clusters:globular, stars: binaries, stars: stellar dynamics 
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1 Introduction 



Recent observational searches suggest that the frequency of primordial bina- 
ries in globular clusters may reach 10% (see Hut at al. 1992 for a review). 
The inferred abundance of binaries is sufficiently large that both dynamical 
and coUisional consequences are important. Several different treatments con- 
clude that binaries are effective in halting core collapse, supporting the core 
and/or driving an expansion phase (Goodman & Hut 1989, McMillan et al. 
1990, 1991, Gao et al. 1991, Heggie & Aarseth 1992). The abundance and 
binding energy distribution of the binaries have a direct impact on observable 
characteristics of globular clusters such as the size of the core radius (Ves- 
perini & Chernoff 1994). The inferred abundance of binaries is consistent 
with that needed by stellar reaction pathways to account for the population 
of millisecond pulsars in the cluster system (Sigurdsson and Phinney 1994). 

The fraction of binaries observed among field stars (> 70% see, e.g., Abt 
1983, Duqucnnoy & Mayor 1991, and Rcipurth & Zinnccker 1993 for the 
frequency of binary stars among prc-main sequence stars) is larger than the 
fraction, hitherto detected, in clusters. The size of the difference may be 
partially due to observational biases: binaries should segregate to the center 
of the cluster and may be difficult to find on account of stellar crowding. 
Or fewer binaries may have formed originally, some may be depleted by 
coUisional effects and some visible members may be swapped by exchange 
reactions for heavy, degenerate and invisible remnants. Many lines of inquiry 
(reviewed by Hut et al. 1992) are being pursued to characterize the present 
binary content and binary binding energy distribution within clusters. 

The key question we will address is whether it is possible for binaries, 
assumed to be primordial, to survive the birth of the cluster. The route 
for the formation of a bound cluster of stars is not well-understood but it 
is plausible that a hydrostatically supported gas cloud is a direct precursor 
to the stellar cluster. Such a cloud might form by coohng instability in the 
Galactic halo (Fall & Rees 1985) or in the cooling region behind a shock 
front after two gas clouds coUide (MacLow & Shull 1986, Shapiro & Kang 
1987). The narrowness of the main-sequence in Galactic globular clusters is 
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evidence that the epoch of star formation is brief. If star formation proceeds 
rapidly compared to the gravitational free-fall time, the resulting collection 
of stellar objects forms far out of virial equilibrium and a period of violent 
relaxation ensues. Stars exchange energy through fast, large-scale variation 
of the gravitational potential. Although the collapse is dissipationless, the 
gross properties (radius, velocity dispersion) of the system undergo damped 
oscillations until a virialized state is reached. During the intervals of high 
stellar density and high velocity dispersion coUisional interactions occur with 
increased rate. 

In this paper we will investigate the fate of a population of primordial 
binaries, paying particular attention to the process of disruption of binaries. 
Since the total binary binding energy can easily exceed the smooth, large- 
scale gravitational potential energy of the stars of the cluster, binary heating 
and softening are also potentially important processes to address. Although 
we treat only point-like stellar objects and assume all gas is removed at the 
conclusion of star formation, our conclusions have some applicability to the 
frequency of collisional interactions of stars with gaseous disks. Elsewhere, 
we will consider the implications for a scenario in which stellar disks are 
common once violent relaxation begins. 

To address the issues, we explore the rate of collisional interactions be- 
tween single stars and binaries. We determine how the binary distribution 
function is altered including the fraction of binaries destroyed and identify 
a characteristic binding energy below which the initial distribution is trun- 
cated. To anticipate our results to some extent, we find that collisional 
effects during the epoch of violent relaxation are small in the following sense: 
in the final, virialized state only soft binaries have been strongly effected (i.e. 
heated or ionized) when the number of stars is as large as it typically is 
for a globular cluster. We explore the scaling of this conclusion with A^. 

The scheme of the paper is the following: in §2 we estimate the binary 
destruction that occurs during the violent relaxation; we discuss the rationale 
for carrying out an N-body experiment; in §3 we describe several tests of the 
N-body code designed to verify that it performs adequately in all regimes of 
interest; in §4 we describe the N-body runs and analyze the results; in §5 we 
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report a refined analytic estimate of tlie cliange in tlie binary distribution 
function. We are able to explain most of the significant trends in the N-body 
simulations and to verify the key role that ionization plays in the numerical 
experiment. In §6 we summarize the main conclusions and describe future 
work. 



2 Violent relaxation and binary destruction: 
a simplified model 

Let ^ = T/|y I be the initial ratio of kinetic to potential energy of the system. 
For the estimates in this section, we represent a stellar cluster by a uniform, 
spherically symmetric density distribution of radius R. If < 0.5 the system 
collapses and oscillates with decreasing amplitude as it virializcs. The first 
cycle of collapse and rc-cxpansion ("bounce") produces the densest central 
conditions. The bounce occurs when the kinetic energy associated with the 
particles' dispersion becomes comparable to the potential energy. iV-body 
calculations show that the discreteness may play a role in determining the 
cluster's size at maximum collapse Rmc (Aarseth et al. 1988). If the initial 
velocity dispersion exceeds a given amount {^^N'^/^) then the minimum size 
scales like Rmc ^ ^Ri- If the system is very cold, so that ^^N~^^^, and if the 
initial positions arc randomly realized for a uniform density distribution then 
the discrete fluctuations induce a background of perturbations and Rmc ~ 
N-'/'Ri. 

If no star gains enough energy to escape and if the cluster remains homo- 
geneous, then the virializcd state satisfies energy conservation so that 

= (1) 

5 Ri ^ ^' 10 Rf ^ ^ 

where M and Ri are respectively the mass and the initial radius of the system. 
Hence, the final radius is 
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If some of the stars escape from the system, changing the mass by AM = 
Mf — Mi and the energy by AE = Ef — Ei (where final values are indicated by 
"f") but leaving the cluster homogeneous, then energy conservation implies 



The final radius is 



10 Rf 5 Ri 

^ 2(1-0 (l + AE/E,) ■ ^ ' 

If the system cools {AE < 0) and/or loses mass (AM < 0) then the final size 
is smaller than the case in which there is no change of energy; if the system 
heats {AE > 0) but does not lose mass, the opposite is true. 

For a system of single stars, van Albada (1982) showed that the magnitude 
of the mass and energy changes during violent relaxation is correlated to how 
cold the system is at the beginning. Small ^ implies violent, large amplitude 
oscillations in the cluster and large |AM| and |A£'|. When binaries are 
present and suffer coUisional interactions within the cluster, the gain or loss 
of their binding energy is accompanied by changes in stellar translational 
energy. 

The amount of energy tied up in the binary population can be substantial. 
Let the number of single stars, binary systems and total stars be Ns, Nh and 
N — Ns+2Nb, respectively. The ratio of the binary internal binding energy to 
large-scale cluster binding energy {SGN^m^/bR) is a = {5fb/6N){R{l/a)), 
where fb — Ni,/N is the fraction of binaries, m is the mass per star and (1/a) 
is the cluster average of the inverse binary semi-major axis. For a tidally 
limited cluster, on a circular orbit of Galactocentric radius Rg with rotation 
velocity Vrot, 

For a newly formed cluster, the most uncertain element in the above estimate 
is probably (1/a). The field binary distribution is peaked about a period 
P ~ 10^'^ d; the period distribution dNb{P)/d\ogP varies hke a Gaussian 
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with (Tiogp = 2.3 (Duquennoy & Mayor 1991). The peak of the field star 
distribution corresponds to a ~ 30 AU but the quantity (1/a) is efi^ectively 
dominated by binaries of the smallest separation. For a cluster born with 
field-hke binaries, taking a to be of order 30 AU and ~ 0.1 (as observed 
today in globulars, rather than ~ 0.6 as observed in the field, Duquennoy 
& Mayor 1991 and Reipurth & Zinnecker 1993) gives an ultra-conservative 
lower limit for a of 0.4. A more plausible estimate is that a extends to ~ 1 AU 
(as observed in the field population, and as required by millisecond pulsar 
production pathways that involve binary mass transfer), implying a ~ 11. 

For the available binding energy to play a dynamical role in the cluster 
evolution, binaries must interact coUisionally. CoUisional interactions occur 
most rapidly at the point of maximum contraction of the cluster when the 
background stellar density and velocity dispersion are largest. In the esti- 
mates here we will focus on ionization, however, our A^-body treatment is 
completely general (§3 and 4) and includes all the point-hke interactions of 
single stars and binaries. We characterize the binary destruction by the "cut- 
off" semi-major axis acut which corresponds to the smallest value of the initial 
semi-major axis a for which the probability of destruction is of order unity. 
Ionization in single encounters and cumulative energy changes by successive 
encounters both contribute to changes in the form of binary binding energy 
distribution and in the cluster's translational energy. 

The change in the number density of binaries (n^) by ionizing collisions 
with single stars of density Ug is 

^ = -nsUbicrv) (6) 
at 

where (crv) is the rate coefficient, averaged over the relative velocity dis- 
tribution. If the single stars and binaries both have Maxwellian velocity 
distributions, the rate coefficient may be written 

(av) = naWthTZix) (7) 

where TZ{x) is a function of the hardness factor x = e/mVg (e is the binary 
binding energy Gm^ /2a and Vg is the single star, one-dimensional velocity 
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dispersion) given by fits to numerical simulations (Hut & Bahcall 1983) as 

" (l + 0.2^/x)(l + exp[x/^])' 

Vth — 3{A/2y/'^Vs is the dispersion in relative velocity and A is a constant. 
For the case of energy equipartition A = 1; for velocity equipartition A ~ 
4/3. The process of violent relaxation leaves the stars in the cluster close to 
velocity equipartition. 

Let the half-mass radius at the point of maximum contraction be Rh,mc- 
The density p ~ MR^^^^, the velocity Vs ^ {GM / Rh^^^f^ and the bi- 
nary hardness x ~ Rh,mc/o-N. The lapse of time spent by the system 
in its maximally contracted form is proportional to the dynamical time 
tdyn,mc ~ {Gpmc)~^^'^- The fractional decrease in binaries of size a is 

S{a) = —^^tdyn,mc (9) 

Ub at 



\^h,mcj 



aN 



(10) 



where C ~ 1.16 and B ~ 1.56 are numerical coefficients obtained by keeping 
all the appropriate factors in the above estimates. If the binary is soft at the 
point of maximum collapse, then we have the following simple expression for 
the fractional decrease 

5 = 5.56-^. (11) 

^h,mc 

A more complete analysis is possible and is given in Appendix A. The cutoff 
semi-major axis (for which 5 ~ 1) is closely tied to the minimum scale of the 
collapse. 

Aarseth et al. (1988) characterized the minimum scale for collapse of 
a cold spherical system without angular momentum. In basic agreement 
with their results our numerical work (in succeeding sections) shows that the 
system collapses to a minimum scale set by root-N fluctuations in the initial 
conditions. The theory predicts Rh,mc oc N~^/'^ and we find 
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For the homogeneous system of interest here, Ri = 2^^^Rh,i, and 



^h,mc l_427V0-36 

which imphes 

R- 

^cut ~ 0.13^^. (13) 

The hardness of this binary, measured in the final viriahzed system, assuming 
no mass or energy loss during the epoch of violent relaxation, and making 
the approximation ^ = is 

9 6 

^-t ~ (14) 

If N> 10"^ only soft binaries are ionized. The simple analysis indicates that 
violent relaxation truncates the binary distribution at a characteristic soft 
scale. 

The estimates based on the simple model are potentially misleading be- 
cause the model has the following shortcomings: 

(1) The model assumes homogeneity while numerical simulations (see, e.g., 
van Albada 1982, Aarseth et al. 1988) show that the collapse of an initially 

homogeneous cluster gives rise to a corc-halo structure when the point of 
maximum contraction is reached and that strong density variations persist 
in the final state. Since the collisional reaction rates depend sensitively on 
density it is not clear how the estimates will be altered. 

(2) The rate coefficient for the model estimate of ionization is based on 
Maxwellian velocity distributions. Since there exist strongly ordered motions 
during the initial phase of violent relaxation, not thermal distributions, the 
rate of binary destruction may be inaccurate. 

(3) In a cold initial state, particles are correlated over long distances; bound 
binaries, triples, etc. are present in the initial conditions. The rate of destruc- 
tion of these structures cannot be described with reaction rate coefficients 
for free particles impinging on a binary as assumed in the model. 

iV-body simulations arc described in §3 and §4 which arc free of these 
shortcomings. Qualitatively, they validate the conclusion of the simple model 
that only soft binaries are ionized. Quantitatively, they give a more accurate 
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description of the processes that are important in the evolution of the binary 
distribution function. 

3 Tests of the numerical code used 

Simulations have been carried out by a special A'"-body code designed for 
coUisional problems (Jayaraman & Chernoff 1992). Forces are calculated us- 
ing an oct tree and a low-order multipole expansion of the particles within a 
cell (following Salmon 1991). Particles have independent time steps; polyno- 
mial approximations are used to specify positions and velocities within given 
timc-intcrvals. Updated polynomials are found using a predict or- corrector in- 
tegration scheme with accuracy 0(At^'^). The tree's structure is exact with 
respect to current particle positions (i.e. it is not allowed to deform as in the 
treatment of McMillan and Aarseth 1993). The accuracy of the integration is 
controlled by two tolerance parameters, one for the integration step size (for 
both relative and absolute error controls) and one for the neighbor sphere 
size. Well-isolated, bound stellar pairs are treated in the following special 
manner. The internal degrees of freedom are regarded as unperturbed Ke- 
pler orbits and the pairs interact with the rest of the system as single entities 
having the total mass and momentum of the constitutents. All non-isolated 
stellar pairs arc treated as individual stars. The classification into isolated 
and non-isolated pairs is updated on each timestep. Close-encounters be- 
tween singles and/or binaries are identified making use of the tree structure. 
When a close-encounter is detected the particles involved are advanced on a 
collision-based timescale, which is typically much shorter than the timestep 
a single star would take in a coUisionless simulation. Since each object is 
advanced independently, the extra cost for handling close encounters in this 
way is not prohibitive. The code maintains coUisionless equilibria including 
Plummer laws. King models and polytropes for many characteristic dynam- 
ical times. 

Several simulations were designed to test critical aspects necessary for the 
simulations of interest. In the first set of tests the collapse of a homogeneous 
sphere of particles was followed. The evolution provides a good illustration 
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of the behavior of the system in the absence of primordial binaries and allows 
calibration of errors and comparison with previous results. 

Table 1 summarizes the results of three runs with N — 1000. The stars 
are distributed homogeneously within a sphere with virial ratios ^ = 0.25, 
0.1 and 0.01. The particle velocities are drawn from a flat 3d distribution 
function whose maximum value is determined by ^. 

The units of all runs have been chosen so that G = 1, the total initial 
mass of the system Mj = 1 and the initial radius of the system Ri — 1. The 
timescale for the homogeneous sphere to collapse is 



The system is already close to the state of virial equilibrium ai t — 1.3; runs 
were stopped at t = 2. Figure 1 provides snapshots of the mass density 
profiles at t = 0.5, 1.1, 1.5, and 2 for the case ^ = 0.01. In agreement with 
van Albada (1982) the snapshots clearly show the development of a core-halo 
structure during the collapse. The central density reaches its maximum value 
at i ~ as expected. 

Table 1 shows, in agreement with van Albada (1982; see his table [1]), 
the physical trends that clusters with smaller ^ have (i) a larger fraction of 
escaping stars, (ii) a greater energy change in the bound stars and (iii) a 
larger ratio of final to initial central densities, Pc,f/Pc,i- In agreement with 
the analytic model of §2 the final to initial half-mass radius, Rhj/Rh,i may 
be predicted from the observed mass and energy changes. Table 1 also gives 
an indication of the size of the numerical errors. For all runs herein, the 
numerical parameters have been adjusted to give relative energy errors ^1% 
over the entire course of the calculation. 

The second test assesses the abihty of the code to track coUisional encoun- 
ters not only during the phase of violent relaxation but also over relaxation- 
length timescales. The single and binary stars are homogeneously distributed 
in a sphere with ^ = 0.01 as above. Here N — 1000 with A^^ — 500 and 
Nh = 250. The binaries have semi-major axis a ~ 9 x 10~^ and zero eccen- 
tricity. In the virialized state, assuming no energy or mass loss, the binary's 
hardness is a; = b/iNa ~ 1.4. 




(15) 
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The evolution of the system was followed well beyond the end of violent 
oscillations until t = 8 (approximately ~ This corresponds to ~ 

WOtdyn: where we employ the definition of Spitzer and Hart (1971) 

_ 1.58i?f 

evaluated in the virialized state. (Specifically, at t — 2, the measured value 
of the Rh implies tdyn ~ 0.07t//.) The half-mass relaxation time trh — 
Ntdyn/26 log(0.4A^) ~ 15tdyn- The long period of quiescent evolution permits 
a significant number of binary-single interactions. 

Figure 2 shows the energy budget of the stars within the systems split 
into four separate categories (analogous plots are shown in Heggie & Aarseth 
1992) (i) into internal energy (i.e. negative binding energy of a binary) versus 
external energy (i.e. translational energy of a single star or of a binary center 
of mass plus large-scale potential energy) contributions and (ii) into retained 
(bound to the cluster) versus escaping particles. The separate quantities 
plotted are labeled Ei^t, Eint,esc, Eext, and Eext,esc- The total energy of the 
cluster is Eint + E^xt] the conserved energy is the sum of all four components. 

Evolution proceeds through several distinct phases. Infall occurs for 
t < tff] there are no escaping particles yet. The binding energy associ- 
ated with the binaries increases slightly (line of medium dashs) because the 
binaries are hard relative to the cold, background stellar distribution. The 
interactions that take place during this phase depend, in part, upon pre- 
existing bound configurations that have the opportunity to collapse and in- 
teract slightly before tff. The coUisional interactions lead to a small degree 
of binary hardening and a concomitant heating of the stellar distribution, 
seen as a slight increase in the external energy (line of short dashes). 

Mass ejection occurs near tff at the point of maximum contraction. This 
is the largest change that occurs over the course of the entire simulation. For 
the binary distribution realized, the external energy dominates the internal 
energy for the escaping matter, but this need not be true in general. The 
figure clearly indicates the corresponding decrease in the external energy of 
the matter that remains bound. 
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After the ejection has occurred, the binaries within the cluster slowly 
evolve. They soften (their internal binding energy decreases) while the exter- 
nal energy of the cluster decreases. There is little discernible mass or energy 
loss during this phase. (Note, no tidal field was imposed on the cluster.) Due 
to the ejection of mass and energy during the violent relaxation, the hardness 
of the binaries is ~ 0.4, somewhat less than the value 1.4 anticipated on the 
basis of the homogeneous model. Since x^l the binaries soften. 

The main check provided by this test is energy conservation: over the 
course of the entire run the relative energy conservation is ~ 10~^. 

4 A^-body simulations 

4.1 Gross Changes 

In this section we describe a set of N-body simulations designed to answer 
the related questions: (1) Is cluster evolution grossly altered by the pres- 
ence of the primordial binaries during the period of violent relaxation? (2) 
How is the binary distribution modified by collisional interactions? (3) Does 
significant mass segregation of the binary population occur? We take the 
initial stellar system to be spherical, homogeneous and cold (^ = 0.01). We 
make these choices in an attempt to maximize the central density and the 
destruction and/or collisional interactions of the binaries. In another paper 
we will discuss the role of inhomogeneities induced by a tidal field in limiting 
the degree of compression. 

The adopted units are G = 1, total initial mass is Mj = 1, initial radius 
Ri = 1. We first generate a set of initial conditions with single stars (mass m) 
and binaries composed of stars (of mass m) and refer to this as the "normal" 
case. For these runs, the initial binary fraction is = 0.05 and all binaries 
begin with circular orbits. The distribution function for the binding energy e 
is logarithmic, i.e. dN — f{e)de oc de/e. The selected range of binary energies 
corresponds approximately to x e [0.13, 1.3] in a homogeneous, virialized 
state, assuming no energy or mass loss and encompasses the critical cutoff 
semi-major axis acut- All calculations are carried out to i = 1.3 at which 
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point the system is close to virial equilibrium. We also consider two variants 
of the normal case. In one variant, we swap a single star of mass 2m for each 
binary ("inert binaries"). In the other, we swap 2 unbound single stars for 
each binary ("all singles"). The suite of runs are summarized in Table 2. 

Three pair- wise comparisons may be made. To assess the full effect of the 
presence or absence of binaries, we compare the normal run to the one with 
all singles. The difference depends on at least two physical effects present in 
the normal runs but absent in the runs with all singles: mass segregation of 
the binaries and collisional encounters between the binaries and the single 
stars. To gauge the effect of mass segregation alone we compare the runs 
with inert binaries to those with all singles. To judge the role of the binary 
internal degrees of freedom alone we compare inert binary runs to normal 
runs. 

Our results shows that the presence of binaries has only a small influence 
on the overall dynamics. Table 3 compares some key quantities at t — 1.3 
for runs. The relative differences in the total external energy of the bound 
particles, the half-mass radius of the bound particles and the mass of the 
bound particles are tabulated. Averages of these quantities are given when 
multiple realizations of the initial conditions were calculated. The first three 
sections of the table show that only small fractional differences have devel- 
oped. Repeated reahzations suggest that the differences are due to the ran- 
dom sampling of the initial conditions; in any case, no systematic differences 
between runs are observed. 

The fourth section of Table 3 provides quantitative information on the 
mass segregation of actual and inert binaries. The final half mass radii of 
singles, inert binaries and real binaries are compared and some systematic 
trends emerge. In all runs with the inert binaries the heavy particles segregate 
compared to the singles. However, in most runs with real binaries the binaries 
are less concentrated than the singles, a consequence of the binary destruction 
and the collisional recoil that occurs most rapidly at the center. Clearly, 
the final spatial distribution of the binaries is as strongly influenced by the 
collisional encounters as by the process of mass segregation. 

Table 4 examines changes in binary binding energy in the normal runs. 
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(Again average quantities are used where multiple realizations were calcu- 
lated.) Although the fractional change in binary binding energy is sizeable, 
the magnitude of that change is small compared to the energy carried off by 
escaping particles. Taken together, the information in Tables 3 and 4 shows 
that the mass and energy losses from ejected particles dominate the changes 
caused by the presence of the binaries. The systematic influence of binaries 
is small compared to the intrinsic variations in escaping particles. 

The maximum degree of contraction is a key factor influencing the rate of 
collisional interaction of the binaries. Table 5 gives the ratio of the initial half- 
mass radius Rh,i to the half-mass radius at the time of maximum contraction 
Rh,mc for the runs with different A^. The best empirical fits to the degree of 
contraction are powerlaws in the total number of particles A^* = A^s + 
Figure 3 displays log Rh^i/ Rh,mc as a function of log A^^, (mean values are used 
when multiple realizations were calculated). The best fit is 

^"'^ -1.13ArO-36, (17) 



Rh,mc 

in good agreement with the analytical and numerical estimate of Aarseth et 
al. (1988) showing that Rh,i/Rh,mc oc N^^^. 

From Tables 3, 4 and 5 we conclude (i) mass and energy loss is domi- 
nated by ejected particles, (ii) the final structure of the bound system is not 
grossly altered by the presence of binaries, (iii) the maximum compression 
during violent relaxation is not greatly altered by the presence of binaries. 
The particular choice of the form of /(e) does not alter these conclusions as 
long as single-star binary interactions are more important than binary-binary 
interactions. For our goal of exploring the binary truncation process during 
violent relaxation, we are free to choose the form of /(e) that best elucidates 
the truncation mechanism, as we do in the next section. We also conclude 
(iv) the spatial distribution of the binaries is sensitive to the physical effects 
of mass segregation and collisional interactions. Consequently, while (i)-(iii) 

^ We use N to indicate the total number of particles, N ~ Ng + 2Ni,, and iV* to stand 
for the "effective" number of particles N^, = Ns + Nf,. When we compare the results of 
section 2 to the A^-body simulations we use the quantity N^, rather than N in the analytic 
formulae. 
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are independent of the particular choice of /(e) that is not true for the final 
spatial distribution of binaries. For example, if /(e) is dominated by high 
binding energy binaries, the system's evolution tends to the case of the inert 
binaries in which unambiguous mass segregation was evident. 

The conclusions should not be significantly altered if ^ is changed. The 
initial conditions are very cold (small ^), in the sense that discreteness rather 
than initial thermal dispersion governs the maximum compression during re- 
laxation. The cold initial conditions make the epoch of relaxation particu- 
larly violent and maximize the coUisional effects. If warmer initial conditions 
were used, particle ejection would decrease, total mass and energy losses 
would diminish and the binary internal energy changes would become rel- 
atively more important. As far as the coUisional processes are concerned, 
warmer initial conditions give an early epoch that more closely resembles the 
quasi-stationary stage of evolution discussed in §3. 

4.2 Ionization and hardening of the binary population 

The changes to the number of binaries and the binding energy distribution 
are the focus of this section. The initial conditions contain a range of binary 
binding energies [e^i„, Cmax]- At later times, a binary is taken to be any pair 
of particles with total binding energy exceeding e^m- 

Figure 4 shows the time evolution of Niy{t)/Nf,{0) for runs with varying 
total particle number N = 1000 (a), 2000 (b), 3000 (c) and 5000 (d). (Mean 
values are plotted for cases with multiple realizations.) Most of the destruc- 
tion evidently occurs close to the instant the minimum radius is reached 
1.1 < t < 1.2. It also of interest that some binary destruction begins before 
that point. The runs with larger N suffer less of this early destruction. 

Figure 5 shows the distribution of binding energies at the beginning (solid 
hues) and at the end of the simulation (t — 1.3, dashed lines) for the different 
values of N. It is quite evident that the softest binaries in the distribution 
suffer the greatest destruction. At the same time, it is also clear that some 
binary hardening occurs; in particular, all the final states include binaries 

with e > e-rnax- 

Although all the binaries included in the above simulations have initially 



16 



circular orbits, none of our essential results are altered by considering a 
distribution of eccentricities. In fact, if ionization is mainly responsible for the 
observed destruction, as we will show below, no alteration is expected because 
the ionization rate depends very weakly on eccentricity (Hut & Bahcall 1983). 
As a test of this point we have carried out two additional simulations with 
N — 1000 (same initial virial ratio and binary fraction as the normal runs) 
but with all binaries having initially eccentric orbits (e = 0.5). In Figure 6 we 
compare the final distribution of binding energy of binaries from these runs 
with that resulting from simulations with binaries having initially circular 
orbits. No significant differences are evident. 

Figures 7 shows the time evolution of the hardness x calculated using 
the numerically determined central velocity dispersion (the maximum value 
anywhere within the system.) The vertical extent of the line segment indi- 
cates the range of the hardness implied by the initial binary binding energy 
distribution. The figure shows that during the period before the phase of 
maximum contraction, the binaries arc hard so, in a thermal distribution, 
one would expect them to harden, not disrupt. As a result, the early binary 
destruction (when hardening is expected) is noteworthy. Moreover, after the 
transition all the binaries are soft and, in a thermal distribution, would tend 
to become softer. As a result, the binary hardening observed at the end of 
the simulation for the largest initial binding energies is noteworthy. 

The data on the hardening prompts a question: how does one explain the 
degree of hardening that is seen? Although the positions of single stars and 
binaries are randomly picked, the cold initial conditions guarantee that cor- 
relations between particles are "frozen in"; the system must evolve for ~ t// 
before the rate of hardening can be calculated with the thermally averaged 
rate coefficients. In fact, the binary hardening inferred using thermally aver- 
aged rate coefficients (Heggie & Hut 1993) and the numerically determined 
density and velocity dispersion history is completely insufficient to explain 
the observed binary hardening before the period of maximum contraction or 
to explain the total binary hardening observed by the end of the simulation. 

It appears that part of the explanation is that there exist bound triples 
present in the cold initial conditions. At t — some binaries are bound 
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to their closest neighbor; the hardening of such bound systems cannot be 
described using the thermally averaged rate coefficients. To show that such 
systems are expected, consider replacing each binary with a single star with 
mass 2m located at the center of mass of the binary. Write the total energy 
of the system composed of the binary and the nearest neighbor: 

1 2 Gmim2 , s 

^triple — ■^l^'V ^ , (18) 

where nil = Mi/N is the mass of the single particle, 1712 = 2mi is the 
mass of the binary, = (2/3)mi is the reduced mass of the system, d = 
{An/Sy^^Ri/N^^^ is the initial average interparticle separation while v'^ ~ 
1.2^GMi/Ri is the initial mean velocity dispersion. We find 

_GMJ/0^_L24\ 

^triple - ^. 1^ ^ ^5/3 ) ^'^) 

which for ^ = 0.01 and ~ 10^, as in our simulations, is negative. The triple 

begins near apocenter; if half the triplc's period is less than tff, pericenter 
will be reached before the point of maximal contraction. If the minimum 
distance between the binary and the single is less than the semi-major axis 
of the binary, a strong interaction and binary hardening is a likely outcome. 
We estimate Uripie, half of orbital period of the single around the binary star, 
to be 

ttriple ~ '^■Uff. (20) 

It is plausible that some binaries have a close interaction before tff and this is 
probably responsible for the hardening that takes place at early times. These 
encounters cannot cause ionization because the total energy of the triples is 
negative. 

As check of the above arguments we have followed the binaries in the 
run with = 3000. The number of binaries initially present in this run is 
Nh = 150; 70 are initially bound to the closest particle to form a triple system. 
The cumulative distribution function of the distance of the distance to this 
closest particle, d^p, is shown in Figure 8a at three different times before 
maximum collapse {t = 0, 0.5 and 1.). As expected, the distance contracts 
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as the system collapses. In Figure 8b all the values of dcp at t = 0.5 and 
t = 1.0 have been rescaled by the factor Rh,i/Rh{t). About 80% of the triples 
collapse more quickly than the system as a whole. This is roughly consistent 
with the cumulative distribution function of ttripie in Figure 8c showing that 
60% of the triples initially present in the system have tfripie < tff- Thus, 
it appears that a sizable fraction of triple configurations have the time to 
collapse, i.e. to reach pericenter, before maximum contraction. 

Figure 8d shows the cumulative distribution function of the ratio of the 
semimajor axis of the binary, a, to the pericenter of the weakly bound element 
of the triple, p. Nine binaries satisfy p < 2a and Uripie < tjf, conditions 
which guarantee a strong interaction. This supports the hypothesis that 
some hardening occurs because of the interaction of singles and binaries 
initially bound in triple systems and that the hardening can occur before the 
bounce. 

It is important to note that a triple is not necessarily destroyed at t ~ 
tff. Hierarchical triples with p<2.75 — 3.5a (Harrington 1972) are inherently 
unstable and, left alone, will eventually interact. The loosely bound outer 
star is expected to harden the binary when the interaction occurs. Also, if 
the triple encounters a single star it may produce an interaction with the 
binary member. It is possible, but unproven, that nearly all the triples will 
suffer strong interactions. 

When triples are abundant the initial configuration may be more impor- 
tant for binary hardening than subsequent thermal single-star binary encoun- 
ters. For a hard binary, the number of single-star binary encounters with a 
minimum distance of separation <a over the course of the period of violent 
relaxation is ~ Nba/ Rh^mc where Nb is the total number of binaries. Since 
Rh,mc ~ R/N^^^ and since the binary binding energies in the numerical runs 
have been scaled so as to keep x in the final virialized state fixed, a oc R/N, 
it follows that the number of binary-single encounters in the N-body calcu- 
lations ~ Nb/N"^^^. Let A^3 be the number of bound triples that ultimately 
give a strong binary interaction. The ratio of triple encounters to single-star 
binary encounters scales like N^N"^/^ / Nb] if the system is so cold that each 
binary is bound in a triple (eqn. and if each triple suffers a strong 
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interaction, then the ratio increases with A^. In fact, the observed change in 
binding energy increases with in the set of simulations. 

In summary, the cold initial conditions include bound triples (and pos- 
sibly higher-order assemblies) that interact and produce hardening in some 
binaries before the point of maximum contraction is reached. Bound triples 
may dominate the hardening after a free-fall time. 

4.3 Cutoff binary binding energy 

We now investigate the quantitative determination of the cutoff binary semi- 
major axis ttcut, or equivalently, the cutoff binding energy Ccut- We first 
introduce an empirical determination based on comparison of the initial and 
final cumulative distribution function (CDF) of the binding energy. Figure 
9 gives the CDFs (initial and final) normalized by the initial total number 
of binaries. We define the cutoff to be the crossover point between the two 
curves. This definition of the cutoff is not the same that we have used in §2 
or §5 so we discuss its plausibility. 

The cross over point identifies the value of the binding energy tcut such 
that the total number of more tightly bound binaries is the same as in the 
initial conditions {AN — Nf{> e) — Ni{> e) = 0). //it were true that 
ionization was the only process responsible for modifying the binary binding 
energy distribution, this definition would not be very useful because some 
destruction can occur at all energies so AA^ < for all e. A more meaningful 
definition might involve a significant relative change {AN/Ni{> e) = const). 

On the other hand, it is clear from the CDFs in Figure 9 that ionization is 
not the only process modifying the properties of the binary star population. 
Significant hardening occurs for large e (evident from the fact that the final 
CDF is much larger than the initial CDF for e > ecut)- This hardening, which 
we have discussed in the previous section, is caused in part by the bound 
triples in the initial configuration. Such encounters do not lead to ionization 
and the rate of hardening encounters does not appear to be very sensitive to 
e. For soft binaries the rate of thermal ionization scales like a (and Figure 
7 shows that x < 1 at the cluster center) so it is clear that ionization is 
most effective at small binding energies. Hence, the numerically determined 
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crossover point occurs at a point where the rate of ionizing encounters exceeds 
the rate of hardening encounters. 

The numerical results are well fit by the following expression 

1.14 GMf 

^cut = (21) 
Since a = GM^/{2Nh), it follows 

= ]VP^^^- (22) 

This numerical result differs in two ways from the theoretical estimate (based 
on the simple model in §2) which gave 

0.13 

"'^"^ - JpM^i- (23) 

The ratio of the two results is ~ 0.33A^°-^^ for N = 1000 the ratio is about 
7, for N — 5000 about 15. That the magnitudes differ by such factors 
is understandable given the variety of approximations of the simple model. 
The difference in scaling is, however, more of a puzzle. 

There appear to be two possible explanations and they are not mutually 
exclusive. First, in the analytic estimate, the asymptotic form for the ion- 
ization rate 7l{x) for a soft binary was used. As the system evolves, a binary 
goes from being hard to soft and the relevant asymptotic description for TZ{x) 
changes. Details may be found in Appendix A. Second, from the discussion 
of the initial correlations of the cold system, it is clear that the frequency 
of occurrence and the binding energy of the triples both depend upon A^. 
Only if A^ ^ 10^ is the nearest star likely to be unbound to the binary. It is 
not easy to perform an a priori estimate of the scaling when such effects are 
present. Given these uncertainties (in addition to those connected with the 
original definition of the cutoff) we also provide a more detailed calculation 
of the change in the full binary distribution function in the next section. 

The 1-dimensional mean square velocity dispersion in the virialized clus- 
ter is 

2 2 GMf 
~ (24) 
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derived by means of the virial theorem and taking the total energy of the 
system equal to \E\ = GMj/ {5Rhj) (see, e.g., Binney & Tremaine 1987) 
where Mj is the final cluster mass and Rhj is the final half-mass radius. The 
hardness of the cutoff binary is 



-cu. = H = 3-75^#. (25) 



and the scaling of Xcut with depends not only on the variation of ecut 
but also on the mass and energy ejection from the system. The values of 
Xcut for the simulations done are [0.445, 0.342, 0.257, 0.154] for = 1000, 
2000, 3000, and 5000, respectively. Assume (i) that Mf and Rhj are well- 
determined. Table 5 shows that the mass and energy loss do not vary greatly 
over the entire set of runs. And, (ii) that the analytic scaling for acut governs 
for A^ > 5000. Then, the best empirical fit for the cutoff hardness is 

^-t = (26) 

Clearly, as A^ increases, Xcut decreases. For all cases investigated, the 
binaries destroyed in the course of the collapse are unimportant in terms of 
the energetics of the cluster and do not modify the dynamics of the collapse 
nor the final state resulting from it. 

5 Analytical estimate of the effects of ioniza- 
tion on the binary stars 

We begin with an estimate of how the binary distribution function evolves if 
thermal ionization were the only operative process and if no mass segregation 
occurred. 

We tabulated the time evolution (from t = to t = 1.3) of number den- 
sity, Usit), and velocity dispersion, Vs(t), of single particles at a set of 11 
Lagrangian mass radii, r^, from the N-body calculations. We assumed the 
binaries to be associated with fixed Lagrangian radii throughout the simu- 
lation. Using the ionization rate (0) we calculate P(e,rM), the probability 
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of destruction for a binary of energy e at tm- Figure 10 shows the function 
Tm) for the runs with = 1000. The largest ionization probabilities 
for each value of binding energy are, as expected, those calculated at the 
innermost Lagrangian radii where Us and Vg reach the highest values. 
The final binding energy distribution may be approximated 

/(e,i = 1.3) = /(e,i = 0)[l-5(6)] (27) 

where 

J Hbd-'r 

Figure 11 shows the initial and the final distribution function for runs with 
= 1000, 2000, 3000 and 5000 (dots give the final binding energy distri- 
bution function for the A^-body simulation, the dashed line gives the initial 
distribution and the solid line gives the above estimate of the final distribu- 
tion). The binary hardening at the largest values of e is evident from the 
points that lie to the right of the initial range of the distribution. Apart 
from this effect, the agreement between the two calculations is satisfactory. 
Thermal ionization is responsible for the large changes that occur at low 
energies. 

Next, we account for the effect of the triple configurations. Figure 12 
shows the initial and final CDF for the run with = 3000 from A^-body 
data (solid lines) and the analytic estimate (dashed line). For small values 
of binding energy (e/m*cr^ < 11) the analytic CDF is determined by thermal 
ionization, as above. The agreement is good for binaries less bound than 
the crossover point. For higher values of the binding energy we assume that 
each binary undergoes a single hardening encounter with a mean change of 
Ae = 0.4e (Heggie 1975). The resultant final CDF agrees well with the 
numerical CDF. 

Finally, we compare the calculation of ecut and the total fraction of bina- 
ries ionized with the results of the N-body simulations (Tables 6 and 7). To 
find Ecut, we've chosen 6 = 0.21 so that the crossover point of the analytic and 
numerical results agree in the run with A^ = 1000. (Any value of 5 of order 1 
is a priori reasonable.) The cutoff energy and the fraction of ionized binaries 
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both follow in the analytic model, once this choice has been made. Values 
are compared with the values obtained directly from N-body simulations in 
§2. The agreement is quite good and suggests that the cutoff energy and 
binary ionization is controlled by thermally averaged ionization. 

6 Conclusions and Future Work 

We have investigated how a population of primordial binaries is modified 

by an epoch of violent relaxation. We have mostly focused on very cold 
initial conditions because they give rise to the maximum degree of contraction 
and, hence, the strongest coUisional effects. We have used a binary fraction 
fb = 0.05, thought to be typical of today's globular clusters. 

With a combination of analytic and N-body calculations we find the fol- 
lowing: 

1. There is no significant change in the gross cluster properties at the end 
of violent relaxation due to the interaction between internal (binary) 
degrees of freedom and translational degrees of freedom. 

2. There is a characteristic binding energy Ccut such that binaries with 
e > €cut are not significantly disrupted, while for e < Ccut, ionization is 
the main destructive process. 

3. The truncation in the energy distribution of the binary population oc- 
curs for hardness parameter x < 1, as calculated in the virialized sys- 
tem. 

4. As the number of particles increases, collisional effects become less 
important and the cutoff binding energy decreases. 

5. The rate of hardening of binaries with e > Scut observed in the A^-body 
simulations is not described by the thermally averaged rate coefficients. 
Cold initial conditions used for the simulations mean some binaries are 
bound to the closest single particle. An interaction between the single 
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particle and the binary is responsible for at least part of the observed 
hardening. 

The scaling of the cutoff energy and the cutoff semi-major axis with the 
number of particles in the system may be related to the size of the system 
at the moment of maximum contraction. Reasonable agreement between 
A^-body simulations and analytic estimates is found. 

Future work may attempt to increase the fraction of binaries, to incor- 
porate a spectrum of masses, to include gas, and to consider inhomogeneous 
initial conditions. 
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Appendix A 



In deriving expression ([TT| ) for 6{a) we have assumed that the cut-off binding 
energy falls in the very soft regime at the moment of maximum collapse of 
the system; the scaling law ( p!3l) for acut that follows from eq. ([TlD is different 
from that obtained from A^-body simulations. An explanation may be due 
in part to the fact that if the cut-off binding energy does not fall in the very 
soft regime at the moment of maximum contraction, as we have assumed 
above, the approximate expression for Tl{x) we have used is incorrect. In 
the following we will show that depending on the hardness factor of the cut- 
off binding energy at the moment of maximum contraction three different 
scaling laws for acut hold. If we define 

d = a/Rh,mc (29) 

we can write 

S{a) = Ca^NTZ {B / [hN]) (30) 

where B and C are numerical constants introduced in the text. 
It is easy to show that for some range of values of aN , 

{5.56a, if aN oo, 

0.47a2iV, if aN ^ B, (31) 

l.Qa^iVe--^/^^, if dN ^0 

Thus it is clear from the above considerations that, in relation to the 
spacing between the solutions of the equation 6{a) = const, used to define 
the cut-off binding energy (or semi-major axis) there exists three different 
regimes: 

1. N oo: 6 enters the linear regime for very small values of a and does 
not depend on so that acut = const. 

2. 'Intermediate' A^: solutions in the regime where 6{a) ~ a^N and thus 
~acut ~ iV-i/2. 

3. Low A^: a steeper decrease with A^ than acut ~ N~^^'^ is expected. 
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The iV-body calculations, described in the text, with ~ 10'^ — 10"^ have 
a scaling for d that is close to the 'intermediate' N regime (a oc N~^'^^ in 
this range). If, in calculating the solution of S{a) = const., we choose the 
value of the constant 0.003) so to get the correct scaling for the range of 
investigated, we obtain the following estimates for acut 

( 0.19Rh,mcN-°-^^, for < 500 

ttcut = I 0.08Rh,mcN-^/'^ for 500 < < 2000, (32) 

[ 5.39 X IQ-'^Rh^rnc, for A^ ^ oo 

For the range of A^ spanned by numerical simulations (950 < < 4750) the 
curve of the solutions of the equation S{a) = const, is well fitted by 

acut = -^Rh,mc for 950 < AT, < 4750 (33) 

In this way we can obtain the observed scaling but there are still two 
difficulties: 1) in order to get the correct scaling in the range of A^ where it is 
observed in A^-body simulations we have been forced to choose a value of the 
constant that is very small and it is not clear why such a small value should 
be connected with the cut-off energy; 2) the numerical constant in eq. (|3^) 
is much smaller than the value found from N-body simulations. The latter 
point is likely to be due, at least in part, to the very approximate estimates 
for the density and velocity dispersion used to calculate S{a) and for this 
reason a more accurate analytical estimate has been carried out in §5 of the 
paper. 
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Figure Captions 



Figure 1 Mass density profiles for the run with N — 1000 and initial virial 
ratio T/\V\ = 0.01 (no binaries) at (a) t = 0.5, (b) t = 1.13, (c) t = 1.3 and 

(d) t = 2. 

Figure 2 Energy budget for the test run with Ng = 500 single particles and 
Nij = 250 binaries. Total external energy of bound particles (short dashed 
line), total external energy of escaping particles (solid line), total (internal) 
binding energy of binaries bound to the system (dashed line) total (inter- 
nal binding) of escaping binaries (long dashed line) (see text for definitions) . 
Note that the internal binding energy of the binaries is taken here with neg- 
ative sign. 

Figure 3 Ratio of the initial half-mass radius to the half-mass radius at the 
time of maximum contraction. The dashed line is the best-fit line and it is 

proportional to N^'^^. 

Figure 4a-d Time evolution of the number of binaries normalized to the ini- 
tial number of binaries. For the values of N (1000, 2000) for which two 
simulations have been carried out the mean value and the semidispersion as 
error bar are plotted, a) = 1000; b) = 2000; c) iV = 3000; d) iV = 5000 
Figure 5a-d Histograms of the binding energy at the beginning of simulation 
(solid line) and at the end of the simulation (dashed line) (a) N — 1000; (b) 
N = 2000; (c) N = 3000; (d) N = 5000. 

Figure 6 Histograms of the binding energy at the end of the simulations with 
N = 1000 and binaries initially with circular orbits (e = 0) (solid line) and 
eccentric orbits (e = 0.5) (dashed line). Dots have been put on the histogram 
corresponding to the runs with binaries with initial eccentricity e = 0.5 to 
distinguish the two curves where they overlap. 

Figure 7a-d Time evolution of the hardness factor x calculated at the center 
of the system (see text for the definition) for the range of values of binding 
energy put in the initial conditions of N-body simulations. The upper and the 
lower point of each bar are the hardness factor of the upper and lower values 
of binding energy initially put in the systems (a) N — 1000; (b) N — 2000; 
(c) N = 3000; (d) N = 5000. 
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Figure 8a-d (a) Cumulative distribution function of the distance, dcp, be- 
tween the closest particle and the binary forming bound triple systems at 
t = (soHd line), at t = 0.5 (short dashed hne) and ai t — 1 (long dashed 
line); (b) as 10a but with the distances dcp at t — 0.5 and t — 1 rescaled by 
a factor Rh^i/ Rh{t); (c) Cumulative distribution function of Uripie (see text); 
(d) Cumulative distribution function of the ratio of the semimajor axes of 
the binaries bound in triple systems to the pericenter distances between the 
binaries and the singles in the triple systems at t — 0. All the data refers to 
the run with = 3000. 

Figure 9a-d Initial (solid lines) and final (dashed lines) cumulative distribu- 
tion functions of the binding energy of binaries in the N-body simulations 
normalized to the initial number of binaries, (a) N — 1000; (b) N — 2000; 
(c) N = 3000; (d) N = 5000. 

Figure 10 Destruction probability P{rM, e) as function of binding energy and 
Lagrangian radius (see text for the analytical derivation) 
Figure lla-d Final binding energy distribution function from iV-body simu- 
lations compared with the final binding energy distribution function derived 
analytically (solid lines). Dashed lines show the initial binding energy distri- 
bution function, (a) N = 1000; (b) N = 2000; (c) = 3000; (d) A^ = 5000. 
Figure 12 Initial and final cumulative distribution function from N-body 
simulation (solid lines) for the run with N — 3000 and the final cumula- 
tive distribution function from analytical estimates (dashed line) calculated 
as described in the text: 1) e/m*(7^ < 11 ionization only; 2) e/m^a^ > 11 
hardening: initial energy shifted by an amount equal to the mean change in 
energy per encounter ~ 0.4e. 
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Table 1 
Summary of test results 





Ei 




AM/Mi 


Pc,f/Pc,i 


Rhj/ Rh,i 


\AEjE,\ 


0.25 


0.44 


0.44 




19.4 


0.45 


5 X 10-* 


0.1 


0.53 


0.61 


0.15 


115.6 


0.26 


3 X 10-* 


0.01 


0.59 


0.92 


0.18 


1227 


0.17 


2 X 10-* 



Table 1: Summary of test results. All runs have Mj = 1, i?j = 1 and 
= 1000. The final time is t = 2, ,^ is the initial virial ratio T/|\/|, Ei 
is the total initial energy, Ef is the total final energy of particles bound to 
the system, AM/M^ is the fractional mass loss, Pcj/ Pc,i is the ratio of final 
to initial central densities, Rhj/Rh,i is the ratio of final to initial half-mass 
radii, and AE/Ei is the fractional change in energy due to computational 
errors. 
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Table 2 



Initial conditions for N-body 



N 


Normal 


# runs 
Inert Binaries 


Singles 






eccentricity 


1000 


2 


2 


2 


5-50 


5-50 





1000 


2 








5-50 


5-50 


0.5 


2000 


2 


2 


2 


5-52 


2.5-26 





3000 


1 


1 


1 


5-60 


1.66-20 





5000 


1 


1 


1 


5-60 


1-12 






Table 2: Initial conditions for the N-body simulations. All runs have = 1, 
Ri — 1 and an initial virial ratio ^ = 0.01. is the number of stars. Normal 
runs include 0.9N singles and 0.05A^ binaries, inert binary runs include 0.9A^ 
singles and 0.05A'" singles (of mass 2m), and single runs include N singles. 
For runs with binaries, the range of binding energies is given e/ma"^ (and 
e/m^a'^), where m is mass per particle and m* = 10~^, and a = \J^^^ is 
proportional to the initial velocity dispersion. 
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Table 3 
Comparisons 



Af 

i V 






Ef,ib 


1000 


0.029 


0.099 


-0.077 


2000 


0.012 


—0.087 


0.091 


dOOO 


n n /I o 


O.Ooi 


— U.04i 


5000 


0.014 


0.139 


-0.145 


N 

1 V 


Rh,b 


^Rh,b-ib 
Rh,b 


^Rh,ib-s 
Rh,ib 


1000 


—0.043 


—0.064 


0.020 


2000 


0.093 


0.zz5 


—0.170 


ouuu 


n 1 1 9 


n n/i 1 


n 1 /i7 


5000 


-0.003 


-0.153 


0.132 


N 


Mf.t 






1000 


-0.004 


-0.019 


0.015 


2000 


0.005 


0.036 


-0.032 


3000 


0.035 


0.001 


0.034 


5000 


-0.014 


-0.012 


-0.002 


N 


(Rh.f(s)-Rhj(ib)) 


(HHj(s)-nuj(b)) 




Rh,fis) 


Rh.fis) 




1000 


0.17 


-0.15 




2000 


0.23 


-0.38 




3000 


0.16 


-0.11 




5000 


0.13 


0.02 





Table 3: The first three sections report fractional differences in the final en- 
ergy (AEf), half-mass radius (ARh) and mass (AMj) for particles remaining 
bound to the cluster at t = 1.3. The parameters of the runs are listed in Table 
2. Each column applies to a pair of runs identical except for their treatment 
of the binary component. The pairs are identified by the subscripts "b", "s", 
"ib" which refer to runs with normal binaries, with all singles and with inert 
binaries respectively. 

The fourth section reports the fractional differences in the half mass radii 
of separate components in the same run. Rhj{s) — Rh f{ib) is the difference 
between the final half-mass radius of singles and inert binaries in the runs 
with inert binaries. Rhj{s) — Rhj{b) is the difference between the final half- 
mass radius of singles and binaries in the runs with normal binaries. 
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Table 4 




Chanf 


les in bindini 


y energy 












1000 


-0.006 


-0.23 


2000 


-0.007 


-0.29 


3000 


0.0 


-0.02 


5000 


0.005 


0.39 



Table 4: A comparison of average relative energy changes for particles 
bound to the cluster at t = 1.3 for runs with normal binaries. AEb is the 
difference between final and initial internal energy (negative binding energy) 
of binaries bound to the cluster, Eggc is the total external energy carried away 
by escaping particles, and EB,i is the total initial binary binding energy. 
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Table 5 





type 




|AAf|/M, 


\^E\/E, 


1000 


Singles 


13.4 


0.20 


0.71 






12.3 


0.19 


0.60 




Inert binaries 


11.0 


0.20 


0.53 






11.6 


0.17 


0.56 




Normal 


13.7 


0.22 


0.79 






13.7 


0.17 


0.63 


2000 


Singles 


13.1 


0.21 


0.68 






14.7 


0.22 


0.93 




Inert binaries 


17.3 


0.23 


0.99 






15.3 


0.25 


0.95 




Normal 


14.8 


0.21 


0.82 






14.7 


0.21 


0.83 


3000 


Singles 


16.7 


0.25 


1.07 




Inert binaries 


17.2 


0.23 


1.00 




Normal 


19.1 


0.23 


1.17 


5000 


Singles 


18.4 


0.24 


1.20 




Inert binaries 


18.4 


0.27 


1.15 




Normal 


23.8 


0.28 


1.51 



Table 5: The degree of compression, mass loss and energy change in the 
N-body simulations. Rh,i/Rh,mc is the maximum ratio of initial to half-mass 
radii over the course of the calculation. AM is the difference between the 
final and the initial mass and AE is the difference between the final and the 
initial external energy of the particles bound to the system. 
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Table 6 



N 


ecwt(N-body)/m*a^ 


ec«t(analytic)/m*(j^ 


1000 


34 


34 


2000 


14 


11.5 


3000 


9 


9.9 


5000 


5 


5.5 



Table 6: A comparison of the cutoff energy observed in the N-body simula- 
tions and the analytic estimate. 
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Table 7 





D(N-body) 


D (analytic) 


1000 


0.48 ±0.05 


0.51 


2000 


0.48 ± 0.05 


0.37 


3000 


0.42 


0.39 


5000 


0.41 


0.40 



Table 7: A comparison of the fraction of disrupted binaries (D) in the 
N-body simulations and the analytic estimate. 
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